Libraries:

# libraries -------------------
library(dplyr)
library(tidyverse)
library(ggplot2)
library(GGally)
library(glmmTMB)
library(TMB)
library(emmeans)
library(DHARMa)
library(lmtest)

Set working directory

knitr::opts_knit$set(root.dir = '/Users/user/Desktop/Data/Regen_RProj/')

Source sapling data:

source("Scripts/SA_Import.R")

Pivot wider to create dataframe where each row is for one plot and has total details for each species group

sa_merge2 <- sapling_merge %>% 
  pivot_wider(names_from = Species_Groups, values_from = Total_Tally)

Import time since data and add it to the PIRI dataset

time_since <- read_csv("CleanData/Treat_Year_Data.csv")

sa_merge3 <- merge(sa_merge2, time_since, by = 'Site')
#log transform time from treatment data
sa_merge3$l.TFT <- log(sa_merge3$Time_from_Treat)

Run the ‘Add_BA’ script and merge with dataset:

source("Scripts/Add_BA.R")

# merge with sapling dataset -------------------
sa_merge4 <- merge(sa_merge3, prism_BA, by = 'Plot_No')

Run ‘Ground_Data.R’ script and add it to sapling dataset:

source("Scripts/Ground_Data.R")

# merge with sapling dataset -------------------
sa_merge5 <- merge(sa_merge4, ground3, by = 'Plot_No')

Source and add veg data:

source("Scripts/Veg_Data.R")

# merge with sapling dataset -------------------
sa.all <- merge(sa_merge5, veg3, by = "Plot_No")

rm(sa_merge3,
   sa_merge4,
   sa_merge5)

Sapling count data is taken in 25m2 plots; basal area is measured in hectares; veg and soil data is taken a 1m2 plots.

Sapling data will be convered into 1m2 plots in order to compare across and reduce the amount of scales of data collection, reducing it to two: 1m2 plots and per hectare observations

sa.all$PIRI.1m <- sa.all$PIRI/25
sa.all$SO.1m <- sa.all$Shrub_Oak/25
sa.all$Other.1m <- sa.all$Other/25

Create log transformed categories for newly added variables, then select for just the desired variables:

sa.all$l.PIRI1 <- log(sa.all$PIRI.1m + 1)
sa.all$l.SO1 <- log(sa.all$SO.1m + 1)
sa.all$l.other1 <- log(sa.all$Other.1m + 1)

sa.all2 <- sa.all %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI1, Shrub_Oak, l.SO1, Other, l.other1, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)


sa.all3 <- sa.all

sa.all3$l.PIRI <- log(sa.all3$PIRI + 1)
sa.all3$l.SO <- log(sa.all3$Shrub_Oak + 1)
sa.all3$l.other <- log(sa.all3$Other + 1)

sa.all3 <- sa.all3 %>% 
  select(Treat_Type, Region, Site, Plot_No, PIRI, l.PIRI, Shrub_Oak, l.SO, Other, l.other, Time_from_Treat, l.TFT, BA_HA, l.BA_HA, PIRI.BA_HA, l.BA_piri, Mineral_Soil, l.Mineral, Litter_Duff, avgLD, avgLD_l, Veg_Total, l.Veg_Total) %>% 
  arrange(Treat_Type)

Select just for numerical vs log and then look at paired plots:

#not log transformed (1m2)
sa.num <- sa.all2 %>% 
  select(PIRI.1m, S0.1m, Other.1m, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(sa.num)
ggpairs(sa.num, aes(color = Treat_Type))


#log transformed (1m2)
sa.numl <- sa.all2 %>% 
  select(l.PIRI1, l.SO1, l.other1, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(sa.numl)
ggpairs(sa.numl, aes(color = Treat_Type))

rm(sa.num,
   sa.numl)

Can see the correlation coefficients for linear (Pearsons) relationships. None of them appear very strong, except for ones that are analogs (avg LD vs mineral soil exposure; ba/ha vs piri ba/ha)

No significant correlations obvious from plots

Above code and naming conventions are the same as SA_GLMM.Rmd script

Looking at pairwise again but with sa.all3 dataset (untransformed sapling counts)

#not log transformed (25m2)
sa3.num <- sa.all3 %>% 
  select(PIRI, Shrub_Oak, Other, Time_from_Treat, BA_HA, PIRI.BA_HA, Mineral_Soil, avgLD, Veg_Total, Treat_Type)

ggpairs(sa3.num)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggpairs(sa3.num, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

#log transformed (25m2)
sa3.numl <- sa.all3 %>% 
  select(l.PIRI, l.SO, l.other, l.TFT, l.BA_HA, l.BA_piri, l.Mineral, avgLD_l, l.Veg_Total, Treat_Type)

ggpairs(sa3.numl)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

ggpairs(sa3.numl, aes(color = Treat_Type))
## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero

## Warning in cor(x, y): the standard deviation is zero
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

rm(sa3.num,
   sa3.numl)

Again no strong correlations

Going to try some scatterplots:

#PIRI & SO
ggplot(sa.all3, aes(x=l.PIRI, y = l.SO))+
  geom_point()+
  geom_smooth(method = lm)+
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# negative for control, positive for fallrx, flat for other

#PIRI and other
ggplot(sa.all3, aes(x=l.PIRI, y = l.other, fill = Treat_Type))+
  geom_point()+
  geom_smooth(method = lm) +
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# negative for control & fall rx, weak negative for springrx and mow, positive for harvest

#PIRI and veg
ggplot(sa.all3, aes(x=l.PIRI, y = l.Veg_Total, fill = Treat_Type))+
  geom_point()+
  geom_smooth(method = lm) +
  facet_grid(~Treat_Type)
## `geom_smooth()` using formula = 'y ~ x'

# positive for harvest, weak negative for springrx, negative for mowrx, fallrx, and control

#relationships varying by treatment type may be a reason why the model are having a hard time in DHARMa

none of these relationships are very strong

Going to start with data transformed to 1m plot and no treatment type

sa.m1 <- glmmTMB(PIRI ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all2,
                 family = poisson)
AIC(sa.m1) #213.7
## [1] 213.6906
sa.m1_sr <- simulateResiduals(sa.m1, n = 1000, plot = T)

testZeroInflation(sa.m1_sr) #zero inflated

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0534, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m1a <- glmmTMB(PIRI ~ l.SO1 + l.other1 + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all2,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m1a) #215.7
## [1] 215.6909
sa.m1a_sr <- simulateResiduals(sa.m1a, n = 1000, plot = T) #passes

testZeroInflation(sa.m1a_sr) #still zero inflated, lets test on data set without 1m transformations - sa.all3

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0539, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
# Untransformed sapling counts
sa.m1b <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m1b) #215.9
## [1] 215.8987
sa.m1b_sr <- simulateResiduals(sa.m1b, n = 1000, plot = T) #fails ks

testZeroInflation(sa.m1b_sr) #still very zero inflated

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.056, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#going to test models with untransformed counts and see if zero inflation issues resolve as variables decrease?





# test piri ba
sa.m2 <- glmmTMB(PIRI ~ l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m2) #217.1
## [1] 217.0877
lrtest(sa.m1b, sa.m2) # p = 0.07
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + 
##     l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)  
## 1  11 -96.949                      
## 2  10 -98.544 -1 3.189    0.07413 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# test ba
sa.m3 <- glmmTMB(PIRI ~ l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m3) #215.1
## [1] 215.1096
#test mineral soil
sa.m4 <- glmmTMB(PIRI ~ l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m4) #213.2
## [1] 213.2266
lrtest(sa.m3, sa.m4) # p = 0.7
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
##   #Df  LogLik Df Chisq Pr(>Chisq)
## 1   9 -98.555                    
## 2   8 -98.613 -1 0.117     0.7323
#test litter depth
sa.m5 <- glmmTMB(PIRI ~ l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m5) #211.2
## [1] 211.249
#test veg cover
sa.m6 <- glmmTMB(PIRI ~ l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m6) #209.4
## [1] 209.3909
lrtest(sa.m5, sa.m6) #p = 0.7
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   7 -98.624                     
## 2   6 -98.695 -1 0.1419     0.7064
#test other
sa.m7 <- glmmTMB(PIRI ~ l.SO + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m7) #253.7
## [1] 253.7722
lrtest(sa.m6, sa.m7) #p less than 0.001
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.SO + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq        Pr(>Chisq)    
## 1   6  -98.695                                
## 2   5 -121.886 -1 46.381 0.000000000009734 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#test shrub oak
sa.m8 <- glmmTMB(PIRI ~ l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m8) #249.3
## [1] 249.3442
lrtest(sa.m6, sa.m8) #p less than 0.001
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq       Pr(>Chisq)    
## 1   6  -98.695                               
## 2   5 -119.672 -1 41.953 0.00000000009348 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sa.m6 looks best, lets test model
sa.m6_sr <- simulateResiduals(sa.m6, n = 1000, plot = TRUE) #doesn't pass

testZeroInflation(sa.m6_sr) #still fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0842, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#try model with treat type and see if it passes?

sa.m9 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m9) #256.3
## [1] 256.386
lrtest(sa.m6, sa.m9) #p less than 0.001
## Likelihood ratio test
## 
## Model 1: PIRI ~ l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq    Pr(>Chisq)    
## 1   6  -98.695                            
## 2  10 -118.193  4 38.995 0.00000006983 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m9_sr <- simulateResiduals(sa.m9, n = 1000, plot = TRUE) #passes

testZeroInflation(sa.m9_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.9994, p-value = 0.97
## alternative hypothesis: two.sided
testResiduals(sa.m9_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031089, p-value = 0.7216
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.2775, p-value = 0.42
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.54
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000662650602409639 ) 
##                                         0.002008032
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.031089, p-value = 0.7216
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.2775, p-value = 0.42
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.54
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000662650602409639 ) 
##                                         0.002008032
#this model looks great. based on modeling below that starts with TT, I'm going to add veg back into this model and run again

sa.m10 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m10) #218.5
## [1] 218.4779
lrtest(sa.m10, sa.m9) #weird, is now important
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq      Pr(>Chisq)    
## 1  11  -98.239                              
## 2  10 -118.193 -1 39.908 0.0000000002662 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.m10_sr <- simulateResiduals(sa.m10, n = 1000, plot = TRUE) #now fails

testResiduals(sa.m10_sr)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.083632, p-value = 0.001886
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0000000000023763, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.86
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.005070281
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.083632, p-value = 0.001886
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.0000000000023763, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 0.86
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.005070281
## sample estimates:
## outlier frequency (expected: 0.00160642570281124 ) 
##                                                  0

Model has much better fit and zero inflation issues are resolved when treatment type variable is included. Time to do elimination again but with TT

sa.m1b <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m1b) #223.4
## [1] 223.4014
# test piri ba
sa.m2 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m2) #223.9
## [1] 223.9656
lrtest(sa.m1b, sa.m2) # p = 0.1
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.BA_piri + l.Mineral + 
##     avgLD_l + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + l.BA_HA + l.Mineral + avgLD_l + 
##     l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  15 -96.701                     
## 2  14 -97.983 -1 2.5642     0.1093
# test ba
sa.m3 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m3) #222.3
## [1] 222.2865
#test mineral soil
sa.m4 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m4) #220.5
## [1] 220.463
lrtest(sa.m3, sa.m4) # p = 0.7
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Mineral + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + avgLD_l + l.Veg_Total + 
##     offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  13 -98.143                     
## 2  12 -98.232 -1 0.1765     0.6744
#test litter depth
sa.m5 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m5) #218.5
## [1] 218.4779
#test veg cover
sa.m6 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m6) #256.4
## [1] 256.386
lrtest(sa.m5, sa.m6) #p less than 0.001
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq      Pr(>Chisq)    
## 1  11  -98.239                              
## 2  10 -118.193 -1 39.908 0.0000000002662 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#test other
sa.m7 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m7) #217.6
## [1] 217.5681
lrtest(sa.m5, sa.m7) #p = 0.3
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + l.Veg_Total + offset(l.TFT) + 
##     (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1 | 
##     Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  11 -98.239                     
## 2  10 -98.784 -1 1.0902     0.2964
#test shrub oak
sa.m8 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.m8) #215.6
## [1] 215.6341
lrtest(sa.m7, sa.m8) #p = .8
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.Veg_Total + offset(l.TFT) + (1 | 
##     Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  10 -98.784                     
## 2   9 -98.817 -1 0.0661     0.7971
#so with treatment type included, veg cover is important, but other species of saplings are not

#check model
sa.m8_sr <- simulateResiduals(sa.m8, n = 1000, plot = TRUE) #fails

testZeroInflation(sa.m8_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.088, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
#try model with so and other saplings?
sa.m5_sr <- simulateResiduals(sa.m5, n = 1000, plot = TRUE)

testZeroInflation(sa.m5_sr) #both still fail ...

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0802, p-value < 0.00000000000000022
## alternative hypothesis: two.sided

I have a very bad headache, so I’m going to stop now. Im so delicate. Anyway, veg seems to messing some stuff up, even if AIC says its important. Take a look at scatter plot/lm of piri and veg for saplings. Then play around with models that have and don’t have it. Could try assinging something other than 1 for the zero inflation…

Starting analysis again

sa.m10 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.m10) #221.7
## [1] 221.6967
sa.m10_sr <- simulateResiduals(sa.m10, n = 1000, plot = TRUE)

testZeroInflation(sa.m10_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0632, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m11 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Treat_Type,
                 family = poisson)
AIC(sa.m11) #223.6
## [1] 223.6341
sa.m11_sr <- simulateResiduals(sa.m11, n = 1000, plot=T)

testZeroInflation(sa.m11_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0874, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m12 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.m12) #213.9 w/1, 250.5 with region
## [1] 250.5171
sa.m12_sr <- simulateResiduals(sa.m12, n = 1000, plot = TRUE) #fails with 1 as ZI; passes with region as zi

testZeroInflation(sa.m12_sr) #fails with 1 as ZI; passes with region as zi

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99973, p-value = 1
## alternative hypothesis: two.sided

Ok, this is the best fit model so far

sa.m12 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.m12) #213.9 w/1, 250.5 with region
## [1] 250.5171
#I'd like to try adding so, other, and veg back into this model, checking AIC and model fit

sa.m13 <- glmmTMB(PIRI ~ Treat_Type + l.Veg_Total + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.m13) #221.7
## [1] 221.6967
sa.m13_sr <- simulateResiduals(sa.m13, n = 1000, plot = TRUE) #fails

testZeroInflation(sa.m13_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.0632, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
sa.m14 <- glmmTMB(PIRI ~ Treat_Type + l.SO + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.m14)
## [1] 250.9238
lrtest(sa.m12, sa.m14) #p = 0.2
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.SO + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  12 -113.26                     
## 2  13 -112.46  1 1.5933     0.2069
sa.m15 <- glmmTMB(PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
#won't converge w/ region
AIC(sa.m15)
## [1] 255.3798
sa.m15_sr <- simulateResiduals(sa.m15, n = 1000, plot = TRUE)

testZeroInflation(sa.m15_sr)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99931, p-value = 0.988
## alternative hypothesis: two.sided

Models that pass DHARMa tests:

sa.f1 <- glmmTMB(PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.f1) #256.4
## [1] 256.386
sa.f2 <- glmmTMB(PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
#won't converge w/ region
AIC(sa.f2) #255.4
## [1] 255.3798
sa.f2_sr <- simulateResiduals(sa.f2, n = 1000, plot = TRUE) #passes

testZeroInflation(sa.f2_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99931, p-value = 0.988
## alternative hypothesis: two.sided
testResiduals(sa.f2_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.054361, p-value = 0.1054
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0628, p-value = 0.588
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.52
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000642570281124498 ) 
##                                         0.002008032
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.054361, p-value = 0.1054
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1.0628, p-value = 0.588
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 1, observations = 498, p-value = 0.52
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.000642570281124498 ) 
##                                         0.002008032
lrtest(sa.f1, sa.f2) # p = 0.3 - can drop SO
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.SO + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1  10 -118.19                     
## 2   9 -118.69 -1 0.9938     0.3188
sa.f3 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~1,
                 family = poisson)
AIC(sa.f3) #213.9 
## [1] 213.9437
sa.f3_sr <- simulateResiduals(sa.f3, n = 1000, plot = TRUE) #fails

testZeroInflation(sa.f3_sr) #fails

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1.09, p-value < 0.00000000000000022
## alternative hypothesis: two.sided
lrtest(sa.f2, sa.f3) #p is less than 0.001, but the AIC falls so much, weird
## Likelihood ratio test
## 
## Model 1: PIRI ~ Treat_Type + l.other + offset(l.TFT) + (1 | Site/Plot_No)
## Model 2: PIRI ~ Treat_Type + offset(l.TFT) + (1 | Site/Plot_No)
##   #Df   LogLik Df  Chisq     Pr(>Chisq)    
## 1   9 -118.690                             
## 2   8  -98.972 -1 39.436 0.000000000339 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
sa.f4 <- glmmTMB(PIRI ~ Treat_Type + offset(l.TFT) + (1|Site/Plot_No),
                 data = sa.all3,
                 ziformula = ~Region,
                 family = poisson)
AIC(sa.f4) #250.5 with region
## [1] 250.5171
sa.f4_sr <- simulateResiduals(sa.f4, n = 1000, plot = TRUE) #passes

testZeroInflation(sa.f4_sr) #passes

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 0.99973, p-value = 1
## alternative hypothesis: two.sided
testResiduals(sa.f4_sr) #passes

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.032779, p-value = 0.6584
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.88128, p-value = 0.89
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00100401606425703 ) 
##                                                  0
## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.032779, p-value = 0.6584
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.88128, p-value = 0.89
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 498, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.000000000 0.004016064
## sample estimates:
## outlier frequency (expected: 0.00100401606425703 ) 
##                                                  0
AIC(sa.f4, sa.f2)
##       df      AIC
## sa.f4 12 250.5171
## sa.f2  9 255.3798
#aic is lower for sa.f4 & more df - therefore better model? maybe graph both to see what is going on?

Plot these babies

library(sjPlot)
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
library(sjlabelled)
## 
## Attaching package: 'sjlabelled'
## The following object is masked from 'package:forcats':
## 
##     as_factor
## The following object is masked from 'package:ggplot2':
## 
##     as_label
## The following object is masked from 'package:dplyr':
## 
##     as_label
library(sjmisc)
## 
## Attaching package: 'sjmisc'
## The following object is masked from 'package:purrr':
## 
##     is_empty
## The following object is masked from 'package:tidyr':
## 
##     replace_na
## The following object is masked from 'package:tibble':
## 
##     add_case
set_theme(base = theme_classic(),
          theme.font = 'serif',
          axis.title.size = 1.5,
          axis.textsize.x = 1.5,
          axis.textsize.y = 1.5,
          title.size = 2.5,
          title.align = "center",
          legend.pos = "right",
          legend.size = 1.5,
          legend.title.size = 1.5,
          #legend.bordercol = "black",
          legend.item.size = .75)

plot_model(sa.f4, 
           type = "pred",
           terms = "Treat_Type")
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

plot_model(sa.f4, 
           type = "pred",
           terms = c("l.TFT", "Treat_Type"),
           axis.title = c("Time from Last Treatment (log transformed)", "Total Count of Pitch Pine"),
           title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           show.zeroinf = T,
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

# plot model two
plot_model(sa.f2, 
           type = "pred",
           terms = c("l.TFT", "Treat_Type"),
           axis.title = c("Time from Last Treatment (log transformed)", "Total Count of Pitch Pine"),
           title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           show.zeroinf = T,
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

plot_model(sa.f2, 
           type = "pred",
           terms = c("l.other", "Treat_Type"),
           axis.title = c("Other Saplings (log transformed)", "Total Count of Pitch Pine"),
           title = "Predicted Counts of Pitch Pine Saplings \n >/= 2.5 cm & < 10 cm DBH",
           legend.title = "Treatment Type",
           line.size = 1,
           show.zeroinf = T,
           ci.lvl = NA,
           colors = c("#D8B70A", "#02401B", "#A2A475", "#81A88D", "#972D15"))
## Warning in checkTerms(data.tmb1$terms, data.tmb0$terms): Predicting new random effect levels for terms: 1 | Plot_No:Site
## Disable this warning with 'allow.new.levels=TRUE'

Generally I’ve been starting with constructing full models that don’t include treatment type. Maybe look at by region as well? Think about random and fixed effects

sa_lr1 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
                  data = sa_lrdata,
                  family = poisson)
AIC(sa_lr1) #209.5

sa1_sr <- simulateResiduals(sa_lr1, n = 1000, plot = TRUE) #does not pass
testZeroInflation(sa1_sr) #yes

sa_lr2 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
                  data = sa_lrdata,
                  ziformula = ~1,
                  family = poisson)
AIC(sa_lr2) #211.5

sa2_sr <- simulateResiduals(sa_lr2, n = 1000, plot = TRUE) #fails
testZeroInflation(sa2_sr) #yes

sa_lr3 <- glmmTMB(Total_Tally~ l.BA_HA + l.BA_piri + l.Mineral + avgLD_l + l.Mineral + (1|Site/Plot_No) + offset(l.TFT),
                  data = sa_lrdata,
                  ziformula = ~Region,
                  family = nbinom2)
AIC(sa_lr3) #215.3

sa3_sr <- simulateResiduals(sa_lr3, n = 1000, plot = TRUE) # fails
testZeroInflation(sa3_sr) #yes

#for zero inflation, tried ~1, Region, Treat Type, l.TFT, 1+Region, 
#nbinom won't work

#maybe I should eliminate some elements from the model and take another look via DHARMa and see if that makes a difference